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ABSTRACT 


The  primary  objective  of  this  work  is  to  improve  the  accuracy  of  discrimination  of  nuclear  explosions  from 
earthquakes,  by  developing  and  implementing  a  transportable  path  correction  algorithm  for  the  Lg/Pg 
discriminant.  A  secondary  result  is  the  development  of  a  depth  discriminant  based  on  Eg  group  velocity. 

A  key  to  understanding  Eg  and  Pg  amplitude  changes,  which  can  be  dramatic  over  just  10s  of  kilometers,  is 
using  observations  of  such  changes  over  short  distances.  We  analyze  changes  in  Eg  and  Pg  along  short 
paths  between  southern  California  seismic  network  (SCSN)  stations  for  relationships  between  propagation 
and  statistics  of  topography,  gravity,  and  crustal  thickness.  A  second  key  element  of  our  approach  is  the  use 
of  categorical  statistical  techniques,  which  frees  us  from  the  common  and  generally  unsupported 
assumption  of  a  linear  relationship  between  geophysical  parameters  and  propagation.  The  principle 
underlying  this  approach  is  that  distinct  crustal  or  path  types  that  have  distinct  effects  on  propagation  can 
be  characterized  by  geophysical  parameters.  The  path  types  can  include  directionally  dependent 
parameters.  Corrections  based  on  calibrated  path  types  are  then  transported  to  uncalibrated  areas  with 
similar  path  types.  This  approach  complements  empirical  methods  that  reduce  the  variance  of  existing  data. 

We  use  an  inversion  approach  to  extend  the  path  correction  techniques  developed  in  southern  California,  to 
global  data  reported  in  the  Reviewed  Event  Bulletins  (REB),  while  maintaining  the  key  elements  of 
identifying  path  types  and  using  short  path  segments.  We  divide  long  source-station  paths  into  short 
segments,  with  each  segment  identified  with  a  particular  path  type,  and  then  invert  for  an  Lg/Pg  ratio 
correction  factor  for  each  path  type.  This  performs  comparably  to  other  waveguide  method  approaches  for 
single  stations,  but  has  the  advantage  of  being  applicable  to  multiple  stations  over  a  broad  region. 

Records  with  Eg  or  Pg  below  the  background  noise  level  contain  important  information  regarding  phase 
blockage.  We  use  such  censored  data  to  determine  path  corrections  and  identify  blockage.  Estimation  of  Eg 
and  Pg  signal-to-noise  (S/N)  levels  based  on  pre-event  noise  is  common,  but  leads  to  inappropriate  use  of 
low  S/N  level  measurements.  The  appropriate  noise  level  to  use  is  that  of  the  pre-phase  noise.  Using  the 
pre-event  noise  underestimates  the  real  noise  level  for  Eg  and  Pg,  with  the  result  that  noise  level 
measurements  of  those  phases  are  used  as  if  they  were  signal.  We  extend  the  inversion  for  path  corrections 
to  use  low  S/N  measurements  as  upper  bounds  on  the  phase  amplitudes  in  a  Maximum  Likelihood  estimate, 
and  apply  it  to  the  global  data  reported  in  the  REB  with  the  result  that  we  clearly  identify  blockage. 

We  also  test  and  model  a  depth  discriminant  based  on  Eg  group  velocity.  We  found  in  southern  California, 
that  Eg  from  earthquakes  and  explosions  differ  in  their  response  to  waveguide  properties.  This  suggests  that 
their  Eg  modal  structures  differ,  which  implies  that  Eg  group  velocities  should  differ  for  earthquakes  and 
explosions.  This  is  confirmed  by  observations  from  shallow  and  non-shallow  sources  in  both  SCSN  and 
Israeli  network  data.  The  SCSN  data  consist  of  582  high  signal-to-noise  (S/N)  Eg  recordings  from  10 
nuclear  explosions,  3  non-shallow  earthquakes,  and  one  very  shallow  earthquake.  The  Israeli  data  consist  of 
630  recordings  from  19  quarry  blasts  and  28  earthquakes,  with  near-surface  to  22  km  hypocentral  depths. 
The  depth  dependent  difference  in  Eg  arrival  times  is  significant  and  so  provides  a  potentially  useful  depth 
discriminant.  We  use  modal  and  finite  difference  synthetics  to  model  the  observations  and  gain  further 
insight  into  their  origin.  Specifically,  we  examine  the  frequency-dependent  effects  on  the  arrival  times  of 
Eg  of  both  modal  excitation  and  Rg-to-Lg  scattering  in  heterogeneous  media. 
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OBJECTIVE 


The  primary  objective  of  this  research  program  is  to  improve  the  accuracy  of  discrimination  of  nuclear 
explosions  from  earthquakes  by  developing  a  transportable  path  correction  algorithm  for  the  Lg/Pg 
discriminant.  A  secondary  result  is  the  development  of  a  depth  discriminant  based  on  Lg  group  velocity. 


RESEARCH  ACCOMPITSHED 


Introduction 

In  this  paper,  we  briefly  review  the  development  and  application  of  path  correction  methods  based  on 
categorical  techniques  and  southern  California  data.  We  focus  on  the  underlying  principles  and  basic 
results,  and  refer  the  reader  to  Baker  and  Eneva  (1999)  and  Baker  (1998)  for  details.  We  then  describe  the 
inversion  approach  used  to  generalize  the  method  for  use  with  long  paths  and  apply  it  to  data  reported  in 
the  REB.  We  discuss  the  importance  of  recognizing  low  S/N  data  by  using  pre-phase  rather  than  pre-event 
noise  levels,  and  the  incorporatation  of  censored  data  into  the  inversion.  We  demonstrate  how  use  of  low 
signal-to-noise  level  (S/N)  measurements  in  a  maximum  likelihood  estimate  identifies  Lg  blockage,  while 
inappropriate  use  of  the  measurements  leads  to  spurious  results,  and  we  present  results  of  inversions  at  a 
single  station  and  over  all  of  Asia.  We  conclude  with  a  discussion  of  a  depth  discriminant  based  on  Lg 
travel  times.  We  present  observations  of  late  arriving  Lg  from  shallow  events  in  southern  California  and 
Israel,  and  mode  and  finite-difference  modeling  of  the  observations  to  assess  different  mechanisms. 

Short  paths  hetween  southern  California  stations,  geophysical  parameters,  and  statistical  techniques 

We  compiled  a  data  set  of  changes  in  Lg,  Pg,  and  the  Lg/Pg  ratio  between  stations  of  the  southern 
California  seismic  network  (SCSN)  together  with  statistics  of  topography,  crustal  thickness,  and  the 
isostatic  gravity  residual.  Statistics  used  included  the  mean,  median,  maximum,  minimum,  roughness 
(various  measures),  average  gradient,  maximum  slope,  and  second  and  third  moments.  Changes  in  absolute 
amplitudes  of  Lg  and  Pg  were  obtained  by  application  of  corrections  for  site  amplifications,  determined  by 
Baker  (1997).  That  permitted  us  to  determine  that  changes  in  Lg  are  significantly  larger  and  are  more 
predictable  than  those  in  Pg  (Baker,  1998). 

We  demonstrated  (Baker,  1998;  Baker  and  Eneva  1999)  that  assuming  the  relationship  between 
propagation  and  various  parameter  values  is  linear  will  cause  some  such  relationships  to  be  missed 
altogether,  while  categorical  techniques  are  appropriate  to  our  assumption  of  the  existence  of  distinct 
crustal  types  (or  path  types  if  we  want  to  use  directionally  dependent  parameters).  We  divided  the  path 
segments  between  stations  into  five  types  using  cluster  analysis  and  the  statistics  of  parameters  listed 
above,  and  found  that  analysis  of  variance  soundly  rejects  the  hypothesis  that  Lg  propagation  differences 
between  the  crustal  types  are  due  to  chance  (Baker,  1998).  Even  simple  divisions  of  the  crust  into  types  by 
a  single  parameter,  such  as  median  elevation  (in  which  case  we  define  low,  medium,  and  high  elevation 
crust),  is  effective  at  reducing  the  variance  in  new  data  collected  in  the  same  area  (Baker  and  Eneva,  1999). 

We  performed  a  test  in  which  we  clustered  short  paths  between  stations  together  with  segments  of  longer 
paths  between  distant  stations.  We  made  corrections  based  solely  on  the  short  paths,  and  applied  them  to 
each  segment  of  the  long  paths.  This  let  us  apply  corrections  based  on  calibrated  areas  (the  short  paths)  to 
some  areas  without  calibration  data  (some  segments  of  the  long  paths)  but  with  similar  crustal  types.  The 
path  corrections  transported  well  from  the  southern  Sierra  to  the  Peninsula  Range,  and  from  other  regions 
to  the  Salton  Trough,  reducing  the  scatter  in  the  Lg/Pg  ratio  in  both  cases.  On  the  other  hand,  corrections 
from  the  southwestern  edge  of  the  Mojave  Block,  when  extended  to  the  central  and  eastern  Mojave  Block, 
actually  increased  the  Lg/Pg  ratio  scatter.  In  that  case,  the  transport  of  corrections  failed.  The  crustal  types 
differ  in  their  effects  on  propagation,  but  the  geophysical  parameters  used  did  not  distinguish  between 
types.  Paleomagnetic  data  indicates  the  far  western  Mojave  Block  has  been  strongly  rotated  relative  to  the 
eastern  block  (e.g.  Dokka,  1997),  and  might  have  been  useful  in  our  study  to  further  distinguish  crustal 
types.  This  points  out  the  need  to  thoroughly  characterize  areas  in  which  corrections  are  to  be  extended. 


An  inversion  extending  use  of  short  segments  of  distinct  path  types  to  long  paths  and  global  data 


Previous  work  on  waveguide  method-based  path  corrections  average  parameters  over  entire  paths  (e.g. 
Rodgers  et  ah,  1999;  Lay  et  al,  1999).  Average  parameter  values  may  however  be  misleading.  For  example, 
a  path  half  in  mountains  and  half  in  a  basin  may  have  the  same  average  elevation  as  a  path  entirely  within 
an  high  plateau,  although  the  3  crustal  types  mentioned  likely  affect  propagation  very  differently.  Thus  it  is 
improtant  that  we  continue  to  base  corrections  on  short  path  segments.  Previous  work  also  has  focused  on 
corrections  at  single  stations.  For  corrections  to  have  a  plausible  physical  basis,  they  must  provide  the  same 
correction  value  for  similar  paths,  even  to  different  stations.  To  retain  the  use  of  short  path  segments  and 
enable  corrections  to  be  made  over  an  entire  region  utilizing  multiple  stations,  we  perform  an  inversion  for 
correction  terms  for  each  crustal  type  based  on  the  equation  ai  =  Cj;*  dj;,  where  ai  is  the  i*  logio(Lg/Pg) 
ratio,  Cj;  is  the  correction  term  for  the  k*  crustal  type,  and  dj^  is  the  distance  of  the  k-  crustal  type  that  makes 
up  the  i-  event-station  path.  This  can  be  written  in  matrix  form  as 
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and  when  there  are  no  censored  data,  can  solved  as  a  simple  least  squares  problem.  We  present  an  example 
of  this  next. 

Figure  1  shows  the  event-station  paths  from  the  June  1996  through  October  1999  REB,  having  S/N  >  2  for 
both  Lg  and  Pg.  For  this  very  simple  test,  we  divide  the  crust  into  types  based  on  increments  of  sediment 
thickness  (<50  m,  50-200  m,  200-500  m,  500-750  m,  750  m-1  km,  1-2  km,  2-4  km,  4-6  km,  6-8  km,  8-10 
km,  and  >  10  km).  Most  of  the  paths  cross  areas  with  less  than  2  km  of  sediments,  and  there  are  virtually  no 
paths  that  cross  sediments  thicker  than  6  km.  so  we  effectively  fit  8  parameters  to  414  Lg/Pg  ratios. 

In  addition  to  the  distance  correction,  we  achieve  an  11%  reduction,  and  a  24%  SMAD^  reduction  (the 
scaled  median  absolute  deviation,  SMAD,  is  the  one-norm  equivalent  of  the  standard  deviation  and  is  a 
more  appropriate  measure  of  scatter  for  these  clearly  non-Gaussian  distributions).  This  demonstrates  that  a 
small  set  of  correction  terms  can  reduce  the  Lg/Pg  scatter  at  many  stations  over  a  large  region. 

Recognizing  and  using  censored  data  to  refine  path  corrections  and  identify  blockage 

When  a  phase  is  blocked,  the  measured  phase  amplitude  will  be  below  the  pre-phase  noise  level.  Such 
records  are  the  most  important  to  identifying  blockage.  When  we  have  such  data,  the  amplitude  measured 
in  the  phase  window  is  used  as  an  upper  bound  on  the  true  phase  amplitude.  Because  we  use  amplitude 
ratios,  we  have  doubly  censored  data.  When  data  are  censored,  we  use  maximum  likelihood  estimation  to 
solve  equation  1.  We  present  two  examples  in  which  this  approach  identifies  areas  of  blockage.  In  the  first, 
we  use  data  at  a  single  station  to  demonstrate  the  errors  that  occur  when  such  censoring  is  not  recognized 
and  dealt  with  appropriately.  In  the  second,  we  demonstrate  the  capability  of  this  approach  to  automatically 
identify  areas  of  blockage  using  multiple  stations  over  a  larger  area. 

There  are  two  important  issues  regarding  censored  data.  The  straightforward  one  is  how  to  analyze  it.  Bias 
due  to  censored  data  has  long  been  recognized  in  seismic  magnitude  estimation  (e.g.  Ringdal,  1977),  and 
methods  of  incorporating  such  data  into  an  inversion  are  well  understood.  The  more  subtle  issue  for  path 
correction  studies  appears  to  be  recognizing  that  data  are  censored.  A  surprising  and  ubiquitous  error  made 
in  this  sort  of  study  is  that  S/N  levels  used  to  select  "good"  data  are  determined  by  comparing  amplitudes 
within  Lg  and/or  Pg  signal  windows  with  pre-event  noise  levels  (e.g.  Xie  and  Mitchell,  1990;  Lay  et  al, 
1999;  Rogers  et  al,  1999).  This  is  an  error.  For  phase  measurements  to  be  plausibly  recognized,  they  must 
have  greater  amplitude  than  the  coda  of  preceding  phases.  Ignoring  this  can  have  serious  ramifications,  as 


the  "Lg"  measured,  for  example,  may  actually  represent  the  amplitude  of  Pn,  Pg,  or  Sn  coda.  This 
undermines  the  basis  for  any  physical  interpretations  of  results  that  may  be  obtained. 


Figure  1:  Event-station  paths  (triangles  denote  stations)  for  events  with  both  Lg  and  Pg  S/N  >  2  as 
reported  in  the  June  1996  to  Oct.  1999  REB,  overlain  on  a  sediment  thickness  contour  map  (contour 
scale  is  in  km). 


In  the  next  simple  example,  we  will  temporarily  neglect  the  issue  of  transportability  in  order  to  make  a 
single  point  clearly.  That  is,  we  address  the  effect  of  using  noise  level  measurements  as  if  they  are  signal 
and  show  how  appropriate  incorporation  of  noise  level  measurements  can  be  used  to  identify  blockage. 
Thus,  although  we  use  the  method  of  dividing  the  crust  into  types,  we  do  so  very  simplistically  for  just  a 
single  station  and  want  to  make  it  clear  that  we  do  not  suggest  that  all  low  elevation  crust  will  block  Lg.  It 
is  just  coincidental  that  for  the  paths  to  this  one  station,  the  areas  that  block  Lg  are  low  elevation. 

A  typical  example  of  waveguide-based  path  corrections,  using  data  from  ZAL,  was  reported  by  Lay  et  al 
(1999),  so  we  perform  a  similar  analysis  using  our  inverse  approach,  this  time  separating  crustal  types  by 
elevation.  We  perform  the  analysis  two  ways.  First,  we  purposely  make  the  mistake  described  above.  That 
is,  we  compare  signal  amplitudes  to  pre-event  rather  than  pre-phase  noise  levels,  so  that  noise  level 
measurements  may  be  used  as  if  they  are  signal  amplitudes.  Then  we  repeat  the  analysis,  but  with  poor  S/N 
records  properly  identified  as  noise  level  measurements. 

From  June  1996  through  October  1999,  the  REB  reported  433  events  with  Pn  S/N  >  2.  Of  those,  only  65 
had  Lg  signal  greater  than  pre-Lg  noise.  Thus,  in  our  first  inversion  85%  of  the  data  input  as  Lg  signal 
measurements  really  repesent  coda  amplitudes  of  earlier  phases.  The  corrections  therefore  are  really  due  to 
fitting  parameters  to  the  ratio  of  Pn  to  Pn,  Pg,  and/or  Sn  coda.  Both  Lay  at  al  (1999)  and  our  first  inversion 
produced  an  approximately  40%  variance  reduction  due  to  the  distance  correction,  but  only  a  few  percent 
more  for  the  waveguide  based  corrections.  Figure  2  demonstrates  why  that  is  the  case.  In  the  top  plot, 
logio(Pn/Lg)  is  binned  for  all  the  data  regardless  of  Lg  S/N  level,  and  plotted  versus  distance  (each  plot 
uses  22  points/bin).  The  large  variance  reduction  with  distance  is  due  to  the  strong  relationship 
demonstrated  by  this  plot.  The  data  however  more  likely  represents  the  ratio  of  Pn  to  Pn  coda,  rather  than 
Pn  to  Lg  amplitude.  The  middle  plot  shows  logio(Pn/Lg)  of  all  433  data  points  again,  plotted  vs.  mean 
elevation.  There  is  no  apparent  relationship,  which  explains  why  the  waveguide-based  corrections  provided 
little  further  variance  reduction.  In  the  bottom  plot,  we  have  separated  the  data  into  two  groups.  Recall  that 
all  the  data  have  Pn  S/N  >  2,  but  now  one  data  set  (the  upper  curve)  contains  the  368  data  points  with  Lg 
S/N  <  2  (where  the  noise  is  pre-Lg  noise).  For  that  data  set  there  is  still  no  relationship  between  mean 
elevation  and  relative  propagation  efficiency.  The  lower  curve  is  from  the  65  data  points  with  Lg  S/N  >  2. 
That  curve  does  show  a  clear  relationship  between  elevation  and  relative  propagation  efficiency.  The 


overall  Pn/Lg  ratio  is  of  course  much  smaller,  as  the  noise  amplitudes  erroneously  used  as  Lg  signal 
amplitude  are  much  larger  than  the  real  signal  amplitudes.  Inversion  results  are  presented  in  Table  1. 


Distance  (km) 


Mean  Elevation  (m) 


Mean  Elevation  (m) 


Figure  2:  Logio(Pn/Lg)  vs.  distance  (top) 
and  logio(Pn/Lg)  vs.  mean  elevation  along 
the  event-station  path  (middle),  both  for  all 
data  regardless  of  Lg  S/N  level  In  these  two 
cases  Lg  amplitude  may  actually  reflect 
amplitude  of  earlier  phase  coda.  The 
bottom  plot  is  logio(Pn/Lg)  for  low  S/N  Lg 
(upper  curve)  and  for  high  S/N  Lg  (lower 
curve).  All  data  are  binned  with  22  data 
points/bin.  The  median  and  two  SMAD 
deviation  confidence  intervals  are  plotted. 
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Table  1:  Results  of  two  inversions  for  path  corrections 
to  ZAL.  In  both,  we  divide  the  crust  into  6  types  based 
on  elevation.  As  described  above  (equation  I),  we 
determine  the  distance  each  path  traverses  each 
crustal  type.  The  important  result  in  this  example  is 
the  correction  term  for  the  lowest  elevation  crustal 
type,  -200  to  300  meters  (top  row). 

In  the  1st  inversion,  we  use  Lg  noise  levels  as  if  they  are 
signal  (second  column).  In  that  case  the  large  number, 
relative  to  other  correction  terms,  (top  row)  erroneously 
predicts  a  small  Pn/Lg  ratio.  In  the  2nd  inversion,  we 
appropriately  use  Lg  noise  levels  as  upper  bounds  on 
signal  amplitudes  in  a  maximum  likelihood  estimate.  As 
an  intermediate  step,  we  determine  the  least  squares 
solution  to  equation  1  for  just  the  high  S/N  data.  There, 
the  correction  term  for  the  lowest  elevation  group  is  zero 
(3'^'^  column),  since  no  paths  with  good  Lg  S/N  traverse 
such  paths.  Finally,  when  we  incorporate  the  censored 
data  into  a  maximum  likelihood  estimate  (last  column), 
the  correction  term  for  the  lowest  elevation  group  is  a 
huge  negative  number,  correctly  identifying  that  crustal 
type  as  blocking  Lg. 

We  are  not  saying  that  all  low  elevation  paths  to  all 
stations,  block  Lg.  Only  in  this  example  using  data 
recorded  at  ZAL  to  demonstrate  the  importance  of 
recognizing  and  appropriately  using  censored  data,  do  all 
the  low  elevation  paths  happen  to  block  Lg. 


The  purpose  of  this  final  example  is  to  demonstrate  that  censored  data  can  be  used  to  systematically 
identify  blockage  within  a  large  region  using  multiple  stations.  To  do  so,  we  use  a  much  simpler  approach 
than  identifying  crustal  types.  We  invert  for  blockage  and  relative  Pn  to  Lg  attenuation  just  for  latitude- 
longitude  grids.  Thus,  for  this  example  we  are  not  addressing  transportability.  This  demonstration  that 
maximum  likelihood  estimation  using  censored  data  can  identify  blockage  over  a  large  region  sets  the  stage 
for  future  work,  which  should  use  more  well-defined  path  types  as  inversion  parameters.  For  this  example, 
our  inversion  parameters  are  simply  geographic  grid  locations,  and  the  data  are  Lg/Pn  amplitude  ratios 
reported  in  the  REB  for  all  the  Asian  stations  from  June  1996  through  October  1999.  We  select  only  data 
where  the  paths  do  not  cross  oceanic  crust.  For  reasons  we  do  not  understand,  the  REB  report  a  large 
number  of  high  S/N  Lg  measurements  from  event-station  paths  crossing  deep  ocean.  We  do  not  address  this 
problem  here,  but  remove  such  data  from  our  inversion  data  set.  This  leaves  us  1834  measurements.  Low 
S/N  level  phase  amplitudes  are  treated  as  upper  bounds  on  the  amplitudes.  We  invert  for  path  effects  for  all 
3°x3°  geographic  grids  that  had  at  least  40  1°  long  path  segments  pass  through  them.  Eigure  3  shows  that 


two  areas  of  blockage  are  identified  (darkest  tiles).  The  most  northeastern  dark  area  covers  a  transition  from 
no  sediment  cover  to  the  east,  to  a  thick  (7  km)  sedimentary  basin  to  the  west.  The  westernmost  dark  area 
corresponds  to  another  thick  (7.65  km)  sedimentary  basin.  Thus,  despite  the  coarse  spatial  averaging  of  this 
simple  example,  two  major  areas  of  blockage  were  automatically  identified. 

^/7’  //7/7’  /^/7’ 


Figure  3:  3x3  degree  grids  through  which  passed  at  least  40,  one  degree  long  path  segments,  along 
with  the  locations  of  1317  events  (dots)  and  18  stations  (triangles)  used  in  the  inversion  The 
important  result  is  that  the  darkest  grid  cells  were  found  to  block  Lg. 


Lg  arrival  times:  a  resurrected  depth  discriminant 

We  noted  earlier  that  Lg  propagation  effects  related  to  waveguide  properties  differ  in  southern  California, 
depending  on  whether  the  source  is  an  earthquake  or  explosion.  The  simplest  explanation  for  this  is  that  the 
shallower  sources  preferentially  excite  shallowly  propagating  modes,  relative  to  the  deeper  sources.  If  so, 
we  should  observe  a  difference  in  Lg  arrival  times.  It  has  long  been  recognized  that  Lg  arrivals  from 
shallow  sources  can  be  delayed  relative  to  Lg  from  deeper  sources.  Pomeroy,  et  al.  (1982),  for  example, 
listed  this  depth  dependence  among  possible  regional  earthquake/explosion  discriminants,  but  abandoned 
the  idea  because  of  the  inconsistency  of  the  observations.  The  data  however  were  analog  and  typically 
around  1  Hz  or  below.  We  find  that  the  discriminant  is  much  more  likely  to  be  successful  today  with 
greater  bandwidth  and  digital  data  permitting  better  waveform  processing  techniques.  Figure  4  shows  two 
waveforms  at  369  km  from  an  earthquake  and  340  km  from  an  underground  nuclear  test  at  the  Nevada  Test 
Site  recorded  at  the  same  station.  The  time  window  defined  by  group  velocities  of  3.6  to  3.0  km/second  is 
marked  on  each  record.  Notice  how  much  later  in  the  window  the  Lg  arrival  from  the  explosion  occurs 
compared  to  the  arrival  of  the  deeper  earthquake. 

Figure  4  also  illustrates  the  method  used  to  measure  Lg  arrival  times.  Unlike  other  seismic  phases,  Lg  is  not 
associated  with  displacements  propagating  along  a  well-defined  ray  path  and  so  there  is  no  distinct  model 
dependent  travel-time  for  the  phase  and  no  travel-time  residual.  Rather,  Lg  is  defined  as  whatever  energy 
arrives  within  a  particular  group  velocity  window,  usually  around  3.6  to  3.0  km/sec,  and  is  modeled  either 
as  shear  waves  trapped  in  the  crustal  waveguide  or  as  higher  mode  surface  waves.  To  quantify  the  relative 
delay  between  seismograms,  we  use  the  median  cumulative  amplitude  (MCA)  of  the  envelope  function  of 
Lg  (Figure  4).  That  is,  the  MCA  is  the  position  that  evenly  divides  the  area  under  the  curve  of  the  envelope 
function  into  two  halves.  Because  the  length  of  the  Lg  group  velocity  window  changes  with  distance,  the 
midpoint  is  normalized  by  the  window  length  and  is  expressed  as  a  fraction  of  total  window  length. 
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Figure  4.  Seismograms  (upper  plot)  from  an 
earthquake  at  369  km  distance  (upper  trace)  and 
a  nuclear  explosion  at  340  km  distance  (lower 
trace),  both  recorded  at  the  same  station,  with 
the  position  of  the  Lg  windows  bracketed  by 
vertical  Ones.  The  envelope  functions  of  just  the 
Lg  windows  are  plotted  versus  position  within 
the  Lg  window  (middle  plot),  with  the  median 
cumulative  amplitude  (MCA)  positions  indicated 
by  vertical  lines.  Better  separation  of  shallow 
events  occurs  when  the  MCA  is  determined  from 
only  those  values  of  the  Lg  envelope  greater  than 
the  mean  value  of  the  pre-Lg  noise  window 
envelope  function  (pre-Lg  window  not  shown). 
The  lower  plot  shows  the  Lg  envelope  functions 
with  that  correction,  and  their  MCAs  (vertical 
lines). 


Earthquakes 
&  Nuclear 
Explosions 
(California) 


1  Shallow 
&  3  Deeper 
Crustal 
Earthquakes 
(California) 


I  T  a 


4-6  Hz  Lg 


200  300  400  500 


Earthquakes 
&  Quarry 
Blasts 
(Israel) 


1-2  Hz  Lg 


200  250  300 


Distiuice  (km) 


4-6  Hz  Lg 


2 


1.5 


^ ^ - 

0  5  10  15  20 

Event  Depth  (km) 

Figure  5.  The  upper  plot  shows  3  examples  in  2 
passbands  in  which  a  measure  of  Lg  arrival 
times  (described  in  Figure  4)  is  plotted  for 
shallow  events  (asterisks)  and  non-shallow 
earthquakes  (circles),  vs.  distance.  The  shallow 
event  Lg  is  always  later.  The  lower  plot  shows  Sg 
minus  Pg  residual  travel  times  vs.  source  depth 
for  133  local  (<  75  km  distance)  Israeli 
earthquakes  and  33  quarry  blasts.  The  median 
earthquake  Sg  minus  Pg  residuals  dencrease 
with  depth  throughout  the  entire  seismogenic 
crust,  relative  to  both  IASPEI91  (used  by  the 
IMS,  squares)  and  the  IsraeU  velocity  model 
(circles).  Quarry  blast  Sg  minus  Pg  residuals 
relative  to  each  model  are  also  shown  (red,  at 
z=0).  This  clearly  shows  that  shear  wave  arrival 
times  depend  strongly  on  source  depth 
throughout  the  crust.  That  is,  the  differences 
observed  at  top  are  not  just  a  near  surface  effect, 
but  rather,  each  increment  of  source  depth 
affects  the  shear  wave  travel  time 


Figure  5  presents  observations  that  demonstrate  the  effect  of  depth  on  S  phase  arrival  times.  We  found  that 
in  all  cases  -  Lg  from  a  very  shallow  earthquake  (1-2  km  source  depth)  vs.  non-shallow  earthquakes 
recorded  at  the  SCSN,  Lg  from  NTS  nuclear  explosions  vs.  non-shallow  earthquakes  recorded  at  the  SCSN, 
and  Lg  from  quarry  blasts  vs.  from  6-22  earthquakes  recorded  in  Israel  -  the  shallow  event  Lg  was  slower. 
No  differences  were  observed  in  Pg  travel  times.  In  addition,  for  near  source  observations  in  Israel, 
earthquake  Sg  was  delayed  relative  to  Pg  as  a  function  of  depth  from  0  to  22  km  deep  sources. 

Mode  synthetic  seismograms  for  sources  from  1  to  30  km  depth,  in  1  km  intervals,  are  shown  in  Figure  6. 
There  is  clearly  more  energy  later  in  the  Lg  window  of  the  shallow  event  seismograms.  This  observation  is 
quantified  in  Figure  7,  where  the  MCAs  of  the  Lg  windows  of  the  synthetics'  envelope  functions  are  plotted 
vs.  source  depth.  There  is  also  considerably  more  energy  in  shallowly  propagating  modes  with  lower  group 
velocities  than  those  that  bracket  the  Lg  window.  In  the  real  Earth,  scattering  at  the  surface  and  by  velocity 
heterogeneities  may  scatter  some  of  that  energy  into  higher  velocity  modes.  That  is  essentially  the 
mechanism  commonly  proposed  for  the  enhancement  of  Lg  by  the  fundamental  mode  (e.g.  Myers  et  al, 
1999).  We  examine  that  mechanism  next  using  finite-difference  synthetics. 

We  use  finite-difference  modeling  to  examine  the  effects  of  scattering  on  Lg.  This  demonstrates  how 
energy  may  be  scattered  from  Rg  to  Lg,  causing  the  observed  delays  of  Lg  for  shallow  events.  We  calculate 
synthetics  for  a  range  of  earthquake  depths  and  distances,  for  both  the  model  of  Figure  6  and  that  same 
model  with  random  heterogeneities  added.  The  heterogeneities  are  zero  mean  Gaussian,  with  a  correlation 
length  of  2.25  km.  and  a  standard  deviation  of  10%  of  the  velocity  in  the  surface  layer,  dropping  to  7%  in 
the  lower  crust.  The  Poisson  ratio  remains  the  same  as  in  each  layer  of  the  homogeneous  model. 

Figure  8  shows  synthetic  seismograms  for  the  homogeneous  and  heterogeneous  models,  for  earthquakes  at 
1  km.  depth  intervals  from  the  surface  to  the  Moho,  at  30  km.  The  Lg  window,  from  3.6  to  3.0  km/sec 
group  velocity,  is  outlined  by  shadowing.  The  homogeneous  model  synthetics  (left)  look  much  like  the 
mode  synthetics,  but  in  this  case  the  fundamental  mode  is  included  (labeled  Rg).  The  seismograms  are 
similar  those  in  Figure  6  (the  mode  synthetics),  but  have  more  detail,  as  the  solution  is  more  complete  than 
that  of  the  mode  summation.  Rg  is  only  generated  by  very  shallow  sources  and  is  much  larger  in  the 
homogeneous  model  case,  as  are  the  higher  modes,  than  for  the  heterogeneous  model  (right).  Seismograms 
from  propagation  through  the  heterogeneous  model  have  the  energy  distributed  more  evenly  throughout  the 
record.  These  simulations  demonstrate  that  Rg  to  Lg  scattering  may  also  be  a  feasible  mechanism  for 
delaying  Lg  arrivals. 

Further  observations  and  modeling  support  the  modal  explanation  for  late  arriving  Lg  from  shallow  events. 
In  particular,  we  observe  better  event  separation  at  higher  frequency.  This  is  predicted  by  the  modal 
explanation,  as  the  modes  are  more  narrowly  constrained  in  depth  at  higher  frequency.  The  Rg  to  Lg 
scattering  model  however,  would  predict  later  Lg  at  lower  frequency,  since  the  higher  frequency  Rg 
scatters  out  away  much  more  rapidly.  Near  source  observations  of  Israeli  quarry  blast  records  confirm  this. 
Below  1  Hz,  Rg  persists  up  to  100  km  distance  from  sources  and  records  tend  to  have  energy  in  either  Rg 
or  Lg,  but  not  in  both.  The  mutual  exclusivity  of  Rg  and  Lg  suggests  that  below  1  Hz,  Rg  may  be  a  major 
contributor  to  Lg.  At  higher  frequencies  however,  where  depth  discrimination  is  stronger,  there  is  no  clear 
Rg  beyond  20  km,  and  no  apparent  mutual  exclusivity  between  Rg  and  Lg. 

CONCLUSIONS  AND  RECOMMENDATIONS 


Studies  using  the  closely  spaced  SCSN  data  identified  two  factors  important  to  successful  transport  of 
waveguide  method-based  path  corrections  for  Lg  to  P  discriminant  values.  One  is  that  parameters  should 
not  be  averaged  over  entire  paths,  but  rather,  long  event-station  paths  should  be  broken  up  into  small 
segments.  The  second  is  that  the  assumption  of  linear  relationships  between  propagation  effects  and 
parameter  values  can  miss  important  relationships,  whereas  categorical  statistics  techniques  allow  us  to 
model  propagation  based  on  the  assumption  that  distinct  crustal  types  exist  and  affect  propagation 
differently  from  each  other.  These  crustal  types  can  be  generalized  to  path  types,  so  that  directional 
parameters,  such  as  whether  the  crust  in  thinning  or  thickening  in  the  ray  direction,  can  be  incorporated. 


An  inverse  method  was  developed  to  retain  these  key  elements  while  determining  path  corrections  for  the 
long  event-station  paths  reported  in  the  REB.  It  was  shown  to  be  effective  in  an  application  to  all  the  Asian 
data  reported  over  3  years.  The  method  reduces  the  Lg/Pg  scatter  for  the  region  using  just  8  parameters 
based  on  sediment  thickness.  More  refined  path  type  separation,  such  as  the  cluster  analyses  performed  for 
SCSN  paths,  should  be  performed.  We  expect  that  would  provide  even  greater  reductions  in  scatter. 

In  the  area  of  path  corrections,  we  recommend  for  all  REB  data,  that  path  types  be  clustered  based  on 
statistics  of  all  possible  geophysical  parameters  and  censored  data  be  identified  and  incorporated  into  the 
maximum  likelihood  inversion  described  here.  This  will  result  in  a  valid  map  of  phase  blockage  most 
relevant  to  IMS  stations.  This  should  be  performed  for  as  narrow  of  depth  range  sources  as  possible,  and 
refined  as  more  accurate  depths  and  more  ground  truth  data  become  available.  Eor  the  path  corrections, 
rigorous  error  analyses  must  be  performed  to  assess  the  accuracy  of  transportability. 

We  have  made  observations  of  delays  in  Eg  and  Sg  from  shallow  events  relative  to  non-shallow  events  in  2 
regions  of  the  world.  The  delays  are  large  and  consistent  enough  to  discriminate  depth,  with  better 
discrimination  at  higher  frequency.  Modeling  indicates  that  the  delay  is  due  to  excitement  of  different 
modes  by  different  depth  sources.  We  recommend  that  an  Eg  group  velocity  discriminant  be  regionalized 
for  all  IMS  stations.  An  effective  depth  discriminant  will  depend  on  accurate  regionalized  travel  time 
models,  and  the  regionalized  discriminant  must  be  rigorously  tested  using  all  available  ground  truth  data. 

Eg  being  composed  of  different  modes  has  profound  implications  for  CTBT  monitoring.  Path  corrections 
for  discriminants,  as  well  as  travel  times,  may  not  be  valid  for  events  at  depths  other  than  those  of  the 
calibration  events.  Ensuring  the  validity  of  any  type  of  Eg  correction  surface  will  depend  on  our 
understanding  of  the  source  depth  effect  on  Eg  and  our  ability  to  more  accurately  assess  source  depth. 


Determination  of  whether  a  measurement  represents  a  signal  or  noise  level  must  be  made  using  pre-phase, 
not  pre-event,  noise  levels.  We  have  demonstrated  that  once  this  is  done,  noise  level  measurements 
incorporated  into  a  maximum  likelihood  estimate  can  automatically  and  accurately  identify  areas  of 
blockage.  We  have  also  demonstrated  how  the  inappropriate  use  of  noise  levels  as  if  they  are  signal 
amplitudes,  a  common  error,  can  easily  lead  to  erroneous  conclusions. 
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Figure  6:  Mode  synthetics  (sans  fundamental  mode)  at  400  km 
distance,  for  a  model  from  the  NTS  to  southern  California  region 
(Asad,  1998)  and  for  0  to  30  km  depth  sources.  Modes  excited  in 
shallow  layers  are  larger  for  shallow  sources  and  so  make  a  greater 
contribution  to  later  Lg.  The  left  plot  shows  seismograms  from  8.0  to 
2.5  km/sec  group  velocities.  The  Lg  window  is  shaded.  On  the  right  is 
just  the  Lg  window,  from  3.6  to  3.0  km/sec  group  velocity. 
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Figure  7:  Median  cumulative 

amplitude  (MCA)  positions  of  Lg 
envelope  functions,  for  the 
seismograms  shown  in  Figure  6). 
Excitation  of  different  modes  is 
consistent  with  observed  delays  of 
shallow-event  Lg. 


WffS  W!^H!C>lffi  ftifi  1 1  i  I  ^ }  ViS 

I  91^  1 1  Cr^RC'i  I 
'j|>B 

^^5  'll  jtJ  i)i|t|T|T^W!ffy[l^j^^ 

^HJ5 1  ^MiS  wj  fi 

^^H[ijr5M  fiHij  jr^Wnfflw^HT^W 

^^^K<S55l 


^^S!wSSv’ 


100  120  140  160  180  100  120  140  160  180 

Time  (seconds)  Time  (seconds) 

Figure  8:  Seismograms  at  400  km  distance  for  the  homogeneous  model  (left),  and  for  sources  from  0 
to  30  km  depth,  from  4.0  to  2.2  km/sec  group  velocity,  with  the  Lg  window  (3.6  to  3.0  km/sec  group 
velocity)  indicated  hy  shadowing.  The  very  large  Rg  for  shallow  events  in  the  homogeneous  model  is 
purposely  clipped,  to  allow  more  details  of  the  smaller  phases  to  he  seen.  The  right  plot  shows  the 
same  for  the  heterogeneous  model.  All  traces  are  plotted  on  the  same  scale. 
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